Group 6 - Final: Woosley Fire

Members: Xinbo Lu, Wenjin Yu, Jody Tran, Doyle Park, Abraham Im

March 13, 2019


Introduction

The Woosley Fire was a recent wildfire that had burned in Los Angeles and Ventura Counties. The fire started on November 8, 2018 at 2:24pm near the Santa Susana Field Laboratory just above the Simi Valley (source). It had burned over 96,949 acres of land damaging structures and forcing evacuations for almost 300,000 people. The cause of this massive fire was due to the Santa Ana winds pushing the fire down south just one the first day. It crossed between San Fernando Valley and the Conejo Valley and then headed towards the Santa Monica Mountains.

Historically, Southern California does not have a history of frequent fires and so recent fires are indicating that these fires are not natural for the ecosystem. Lands that are burned more than once in a twenty year span will allow invasive plants to grow and that causes more wildfire in the future. This effects a lot of the animals that live near the Santa Monica Mountains. Although large animals like mountain lions and deers can escape wildfires, there are smaller animals such as reptiles and amphibians that are unable to do so.

For the people that live in Los Angeles and Ventura Counties, they were forced to leave their homes and were unable to return to their homes until the fire was continuing moving to the west side of the Santa Monica Mountains (source). It was not until November 21, 2018 at 6:11pm that the fire was finally contained.

Purpose of Our Project

Based on different bands, we will find out how can surface water content, vegetation, and elevation can affect the fire. Also, calculate the NBR (normalized burn ratio) for pre fire and post fire, then derive the delta NBR (pre fire NBR – Post Fire NBR) to estimate the fire’s severity. The analysis should improve the evacuation plan, fire hazard prevention, and the reaction plan.

Research Question

With the use of looking at before and after landscapes of wildfire(s), how are we able to utilize different analysis such as vegetation and various indexes to help predict future/potential wildfires?

In [61]:
import os
import rasterio
import pyproj
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import glob
#import earthpy as et
In [8]:
os.chdir('F:\\GEOG458\\Geog458_final\\LC08_L1TP_041036_20181018_20181031_01_T1\\')

Pre Fire Analysis

This section will be looking at the ROI (Region of Interest) in terms of before the Woosley Fire.

Opening the Lidar image

Here we are taking a path of the Lidar image and returning a print of the dataset. To highlight a few parts of the print, the driver is the file format it is taking in which in this case is GTiff, the data image size is 7761 (width) by 7891 (height), and the data type is uint16.

Filtering and Storing

We filter all of the GTiff files and stored them. Here we have two files stored (examples presented) but we will be looking at five out of the twelve GTiff files. These files show the best colors (blue, red, and green) for when we make our visual maps.

In [9]:
## background information about the Lidar image
with rasterio.open('LC08_L1TP_041036_20181018_20181031_01_T1_B1.TIF') as src:
    print(src.profile)
{'driver': 'GTiff', 'dtype': 'uint16', 'nodata': None, 'width': 7761, 'height': 7891, 'count': 1, 'crs': CRS.from_dict(init='epsg:32611'), 'transform': Affine(30.0, 0.0, 263385.0,
       0.0, -30.0, 3948615.0), 'tiled': False, 'interleave': 'band'}
In [71]:
## filter all files end with .tif and store them 
LS8_Bands = glob.glob("F:/GEOG458/Geog458_final/LC08_L1TP_041036_20181018_20181031_01_T1/*.tif")
LS8_Bands[0:2]
Out[71]:
['F:/GEOG458/Geog458_final/LC08_L1TP_041036_20181018_20181031_01_T1\\LC08_L1TP_041036_20181018_20181031_01_T1_B1.TIF',
 'F:/GEOG458/Geog458_final/LC08_L1TP_041036_20181018_20181031_01_T1\\LC08_L1TP_041036_20181018_20181031_01_T1_B2.TIF']

General Overview of Landsat 8 Band 4

Here we are looking at Band 4 image to show the area that we will be looking for our project. This band shows the colors of the pre/post fire really well in grey. But still this does not show much so we will have to create a subset and normalize it.

cite: raster processing

In [74]:
## Prefire Landsat 8 Band 4, full image
## Open Band 4 by using earthpy

import rasterio

with rasterio.open(LS8_Bands[3]) as src:
    LS8_Band4 = src.read(1)

plt.figure(figsize=(15,15))
plt.imshow(LS8_Band4, cmap = 'Greys_r')
plt.title('A general view of Landsat8 - Band4',fontsize=20)
plt.xlabel('Column #')
plt.ylabel('Row #')
Out[74]:
Text(0, 0.5, 'Row #')

Creating a subset of Landsat8 Band 4 ROI (Region of Interest)

Here Band 4 and it is redefined so that the visual is a little better to see the details of the land.

In [106]:
## Create a subset of the full resolution Landsat8 image 
##window = rasterio.windows.Window(1000, 4000, 2500, 2500)
window=rasterio.windows.Window(1500, 5000, 1500, 1500)
with rasterio.open(LS8_Bands[3]) as src:
    subset = src.read(1, window=window)

plt.figure(figsize=(15,15))
plt.imshow(subset, cmap = 'Greys_r')
plt.title('Pre fire Series: Landsat8 Band 4, Subset to ROI', fontsize=20)
plt.xlabel('Column #')
plt.ylabel('Row #')
Out[106]:
Text(0, 0.5, 'Row #')
In [18]:
## create a stack of Lidar Images

landsat_path = "C:/Users/LXB19/Documents/GEOG458/Group6_Final/prefire_output/LS08_PreFire.tif"
land_stack, land_meta = earthpy.spatial.stack(LS8_Bands, out_path= landsat_path)

Stack of twelve bands

Using earthpy, we created a stack of Band 1 through 12. Opening them with rasterio and renaming them so that it will be easier to work with individual bands throughout the project. The color we used is grey so it shows a really good contrast on top of the black background.

Normalizing

We normalize the grid values so that when we make our visual maps, the legend is scaled between 0.0 to 1.0.

Note: we multipled by two so that the colors are enhanced and quality.

In [7]:
## Using earthpy to create a full stack of bands, Band 1 - 12. 

with rasterio.open(landsat_path) as src:
    landsat_prefire = src.read()

band_titles = ["Band 1", "Band 2(Blue)", "Band 3(Green)", " Band 4(Red)", "Band 5(NIR)",
               "Band 6 (SWIR1)", "Band7(SWIR2)", "Band 8(Panchromatic)", "Band 9(Cirrus)", 
               "Band 10(TIRS1)", "Band 11(TIRS2)", "Band 12 (QBA: Quality band)"]

ep.plot_bands(landsat_prefire, title = band_titles, cmap = "Greys_r")
Out[7]:
(<Figure size 864x864 with 12 Axes>,
 array([[<matplotlib.axes._subplots.AxesSubplot object at 0x000002040328C128>,
         <matplotlib.axes._subplots.AxesSubplot object at 0x0000020402F83160>,
         <matplotlib.axes._subplots.AxesSubplot object at 0x0000020402FA96D8>],
        [<matplotlib.axes._subplots.AxesSubplot object at 0x0000020402FD1C50>,
         <matplotlib.axes._subplots.AxesSubplot object at 0x0000020403001208>,
         <matplotlib.axes._subplots.AxesSubplot object at 0x000002040302A780>],
        [<matplotlib.axes._subplots.AxesSubplot object at 0x0000020403053CF8>,
         <matplotlib.axes._subplots.AxesSubplot object at 0x0000020403082278>,
         <matplotlib.axes._subplots.AxesSubplot object at 0x00000204030A9828>],
        [<matplotlib.axes._subplots.AxesSubplot object at 0x00000204030D4DA0>,
         <matplotlib.axes._subplots.AxesSubplot object at 0x0000020403103358>,
         <matplotlib.axes._subplots.AxesSubplot object at 0x000002040312A8D0>]],
       dtype=object))
In [14]:
# Function to normalize the grid values
def normalize(array):
    """Normalizes numpy arrays into scale 0.0 - 1.0"""
    array_min, array_max = array.min(), array.max()
    return ((array - array_min)/(array_max - array_min)) * 2

Pre Fire Visuals: Natural Colors

Here we use Bands 2, 3, and 4 which are blue, green, and red respectively so show the natural colors of before the Woosley Fire.

In [108]:
with rasterio.open(LS8_Bands[3]) as src:
    red_band = src.read(1, window=window)

with rasterio.open(LS8_Bands[2]) as src:
    green_band = src.read(1, window=window)

with rasterio.open(LS8_Bands[1]) as src:
    blue_band = src.read(1, window=window)

Nred = normalize(red_band)
Ngreen = normalize(green_band)
Nblue = normalize(blue_band)

rgb = np.dstack((Nred, Ngreen, Nblue))

plt.figure(figsize=(15,15))
plt.imshow((rgb * 255).astype(np.uint8))
plt.title('Pre Fire Series: Natural Color Image', fontsize=20)
plt.xlabel('Column #')
plt.ylabel('Row #')
Out[108]:
Text(0, 0.5, 'Row #')

Pre Fire: Vegetation

This next visual uses Bands 3 (Green), 4 (Red), and 5 (NIR) to show the vegetation around Malibu before the fire.

In [107]:
with rasterio.open(LS8_Bands[4]) as src:
    red_band = src.read(1, window=window)

with rasterio.open(LS8_Bands[3]) as src:
    green_band = src.read(1, window=window)

with rasterio.open(LS8_Bands[2]) as src:
    blue_band = src.read(1, window=window)

Nred = normalize(red_band)
Ngreen = normalize(green_band)
Nblue = normalize(blue_band)

FlaseColor = np.dstack((Nred, Ngreen, Nblue))

plt.figure(figsize=(15,15))
plt.imshow((FlaseColor * 255).astype(np.uint8))
plt.title('Pre Fire Series: False Color Image to Visualize Vegetation', fontsize=20)
plt.xlabel('Column #')
plt.ylabel('Row #')
Out[107]:
Text(0, 0.5, 'Row #')

Pre Fire: NDVI (Normalized Difference Vegetation Index)

Here is where the normalized function that we did earlier comes to good use. We want to see the NDVI for this area before the fire because NDVI will show how healthy the vegegation was. We use Band 4 (Red) and 5 (NIR) to make this visual. To calculate the NDVI:

 - Band 4 / (Band 5 + Band 4)  

The land before the fire has a range between 0.2 and 0.3 so keep that in mind for when we look at post fire...

cite: calculating NDVI

In [105]:
## calculate NDVI 

with rasterio.open(LS8_Bands[3]) as src :
    prefire_red = src.read(1, window=window)
    

with rasterio.open(LS8_Bands[4]) as src :
    prefire_NIR = src.read(1, window=window)

np.seterr(divide='ignore', invalid='ignore')

ndvi = (prefire_NIR.astype(float) - prefire_red.astype(float)) / (prefire_NIR.astype(float) + prefire_red.astype(float))

plt.figure(figsize=(15,15))
plt.imshow(ndvi,cmap='inferno')
plt.colorbar(shrink=0.5)
plt.title('Pre Fire Series: Normalized Difference Vegetation Index', fontsize=20)
plt.xlabel('Column #')
plt.ylabel('Row #')
Out[105]:
Text(0, 0.5, 'Row #')

Pre Fire: mNDWI (modified Normalization Difference Water Index)

mNDWI is similar in some ways like NDVI but it is primarily utilised for removal of built up land noise. In our case, we wanted to see the water index of the California land before the fire. We use Band 6 (SWIR1) and 3 (Green) to make this visual. To calculate the mNDWI:

- (Band 3 - Band 6) / (Band 3 + Band 6)

We can say that the land already is very dry to begin with.

In [104]:
## mNDWI
with rasterio.open(LS8_Bands[2]) as src :
    green = src.read(1, window=window)
    
with rasterio.open(LS8_Bands[5]) as src :
    SWIR = src.read(1, window=window)

np.seterr(divide='ignore', invalid='ignore')

mNDWI = (green.astype(float) - SWIR.astype(float)) / (green.astype(float) + SWIR.astype(float))

plt.figure(figsize=(15,15))
plt.imshow(mNDWI,cmap='Spectral')
plt.colorbar(shrink=0.5)
plt.title('Pre Fire Series: Modified Normalized Difference Water Index',fontsize=20)
plt.xlabel('Column #')
plt.ylabel('Row #')
Out[104]:
Text(0, 0.5, 'Row #')

Pre Fire: NBR (Normalized Burn Ratio)

Normalized burn ratio highlights burned areas and estimate fire severity also very similar to NDVI. We use Band 5 (NIR) and 6 (SWIR1) to make this visual. To calculate the ratio:

- (Band 5 - Band 6) / (Band 5 + Band 6)

You can basically see how damaged the land is on the visual map because the lower numbers represent heavier burns while the higher means less. The colors are reversed so that we can show the low colors of the burn in the warmer scaler.

In [84]:
## Pre fire, Normalize Burn Ratio 

with rasterio.open(LS8_Bands[4]) as src :
    prefire_NIR = src.read(1, window=window)
    
with rasterio.open(LS8_Bands[5]) as src :
    prefire_SWIR = src.read(1, window=window)
    
np.seterr(divide='ignore', invalid='ignore')

pre_fire_NBR = (prefire_NIR.astype(float) - prefire_SWIR.astype(float)) / (prefire_NIR.astype(float) + prefire_SWIR.astype(float))

plt.figure(figsize=(15,15))
plt.imshow(pre_fire_NBR,cmap='CMRmap_r')
plt.colorbar(shrink=0.5)
plt.title('Pre Fire Series: Normalized Burn Ratio (NBR)', fontsize=20)
plt.xlabel('Column #')
plt.ylabel('Row #')
Out[84]:
Text(0, 0.5, 'Row #')

Post Fire

This section will show the ROI from around January 2018 which is about two months after the fire had occured. Many of the beginning steps such as reading files and normalizing the bands are repeated for the post fire files so let's jump right into the first visual of the post fire!

In [20]:
## user rasterio to open the images for post fire series

os.chdir('F:\\GEOG458\\Geog458_final\\LC08_L1TP_041036_20190122_20190205_01_T1\\')
date = '2019-01-22'

filepath2 = 'LC08_L1TP_041036_20190122_20190205_01_T1_B2.TIF'
with rasterio.open(filepath2) as src:
    blue = src.read(1, window=window)
    
filepath3 = 'LC08_L1TP_041036_20190122_20190205_01_T1_B3.TIF'
with rasterio.open(filepath3) as src:
    grn = src.read(1, window=window)

filepath4 = 'LC08_L1TP_041036_20190122_20190205_01_T1_B4.TIF'
with rasterio.open(filepath4) as src:
    red = src.read(1,window=window)

filepath5 = 'LC08_L1TP_041036_20190122_20190205_01_T1_B5.TIF'
with rasterio.open(filepath5) as src:
    nir = src.read(1, window=window)
    
filepath6 = 'LC08_L1TP_041036_20190122_20190205_01_T1_B6.TIF'
with rasterio.open(filepath6) as src:
    swir = src.read(1, window=window)
In [21]:
## Normalize the bands for natural color image
PostFire_Nred = normalize(red)
PostFire_Ngrn = normalize(grn)
PostFire_Nblue = normalize(blue)
PostFire_Nnir = normalize(nir)

Post Fire: Natural Colors

Here we can see some of the damages near the coastline of the aftermath. But with the natural colors we cannot see how much has been burned away. So we will again be calculating the vegetation, NDVI, mNDWI, NBR, and the dNBR to compare before and after.

In [86]:
PostFire_rgb = np.dstack((PostFire_Nred,PostFire_Ngrn,PostFire_Nblue)) #stack rgb plot

plt.figure(figsize=(15,15))
plt.imshow((PostFire_rgb * 255).astype(np.uint8))
plt.title('Post Fire Series: Natural Color Image',fontsize=20)
plt.xlabel('Column #')
plt.ylabel('Row #')
Out[86]:
Text(0, 0.5, 'Row #')

Post Fire: Vegetation

The next sub-slide will show high damages from the fire in the middle range near the coastline. The small parts that are highlighted as "white" are the cities near the fire. These urban areas where most likely evacuated during the fire.

In [87]:
PostFire_FalseColor = np.dstack((PostFire_Nnir,PostFire_Nred,PostFire_Ngrn)) #stack false color plot

plt.figure(figsize=(15,15))
plt.imshow((PostFire_FalseColor * 255).astype(np.uint8))
plt.title('Post Fire Series: False Color Image to Visualize Vegetation', fontsize=20)
plt.xlabel('Column #')
plt.ylabel('Row #')
Out[87]:
Text(0, 0.5, 'Row #')

Post Fire: NDVI (Normalized Difference Vegetation Index)

Here we will see that the post fire vegetation has extremely worst conditions. The land is basically in the range of 0.2 and this means that the majority of the ROI has very low vegetation.

Note: The tip of the coastline near Malibu is not aligned with each other. Although the code for this is all correct, the plot configuration has slightly shifted it.

In [101]:
def calc(dat1,dat2):
    dat1 = dat1.astype(float)
    dat2 = dat2.astype(float)
    calc_outcome = (dat1 - dat2) / (dat1 + dat2) #ls8 calc equation
    return calc_outcome


ndvi = calc(nir,red) 
plt.figure(figsize=(15,15))
plt.imshow(ndvi, cmap='inferno')
plt.colorbar(shrink=0.5)
plt.title('Post Fire Series: Normalized Difference Vegetation Index', fontsize=20)
plt.xlabel('Column #')
plt.ylabel('Row #')
#plot calculated ndvi
Out[101]:
Text(0, 0.5, 'Row #')

Post Fire: mNDWI (Modified Normalized Difference Water Index)

You will see that after the fire, the heavily middle area of the ROI had a larger difference of difference water index. The areas that are more green means that they had a smaller difference. Remember that earlier in vegetation, the areas that have less vegetation were the urban areas. In here you will see those similar spots will have a smaller water index as well.

In [103]:
mNDWI = calc(grn,swir) #call calc function for mNDWI
plt.figure(figsize=(15,15))
plt.imshow(mNDWI, cmap='RdYlGn')
plt.colorbar(shrink=0.5)
plt.title('Post Fire Series: Modified Normalized Difference Water Index', fontsize = 20)
plt.xlabel('Column #')
plt.ylabel('Row #')
#plot calc mNDWI
Out[103]:
Text(0, 0.5, 'Row #')

Post Fire: NBR & Delta NBR (Normalized Burn Ratio)

Here you will see the major areas of the burn ratio to be low in the NBR visual. We have also done a dNBR as well with the color scale in reserved colors so that the darker purple shade can correlate with the NBR. And the way to calculate the delta NBR is:

- pre fire NBR – post fire NBR

cite: calculating NBR

In [96]:
## Post fire NBR

with rasterio.open(filepath5) as src :
    nir = src.read(1, window=window)
    
with rasterio.open(filepath6) as src :
    swir = src.read(1, window=window)
    
np.seterr(divide='ignore', invalid='ignore')

post_fire_NBR = (nir.astype(float) - swir.astype(float)) / (niv.astype(float) + swir.astype(float))

plt.figure(figsize=(15,15))
plt.imshow(post_fire_NBR,cmap='CMRmap')
plt.colorbar(shrink=0.5)
plt.title('Post Fire Series: Normalized Burn Ratio (Post NBR)', fontsize=20)
plt.xlabel('Column #')
plt.ylabel('Row #')
Out[96]:
Text(0, 0.5, 'Row #')
In [98]:
## Delta NBR 
delta_NBR = (pre_fire_NBR - post_fire_NBR)

plt.figure(figsize=(15,15))
plt.imshow(delta_NBR,cmap='CMRmap_r')
plt.colorbar(shrink=0.5)
plt.title('PostFire Series: Delta Normalized Burn Ratio (dNBR)', fontsize = 20)
plt.xlabel('Column #')
plt.ylabel('Row #')
Out[98]:
Text(0, 0.5, 'Row #')
In [2]:
import numpy as np
from numpy import gradient
from numpy import pi
from numpy import arctan
from numpy import arctan2
from numpy import sin
from numpy import cos
from numpy import sqrt
from numpy import zeros
from numpy import uint8
import matplotlib
import matplotlib.pyplot as plt
import richdem as rd
In [3]:
matplotlib.rcParams['figure.figsize'] = (15, 15)

dem_path = 'DEM_10m/DEM_10m.tif'
dem = rd.LoadGDAL(dem_path)

font = {'family': 'sans-serif',
        'weight': 'bold',
        'size': 16,
        }

Looking at DEM (Digital Elevation Models)


DEM Slope, Aspect, and Hillshade

Digital Elevation Models (DEM) provide a representaion of surface topography (elevation) in 2D space. DEMs are a 3D representaion of a terrain’s surface such as the Earth. For Woosley Fire, we can see that the DEM for this area allows us to make an assumption that low elevation had higher fire damages versus areas that have higher elevation. The areas with higher grounds are the more urban areas.

cite: calculating DEM

In [4]:
plt.imshow(dem, interpolation = 'none')
plt.clim(0, 2500)
plt.colorbar(shrink = 0.5)
plt.title("Digital Elevation Model", fontdict = font)
plt.show()
In [5]:
slope = rd.TerrainAttribute(dem, attrib = 'slope_riserun')
plt.imshow(slope, cmap = 'magma')
plt.clim(0, 95000)
plt.colorbar(shrink = 0.5)
plt.title("DEM Slope", fontdict = font)
plt.show()
In [7]:
aspect = rd.TerrainAttribute(dem, attrib = 'aspect')
plt.imshow(aspect, cmap = 'jet')
plt.clim(10, 360)
plt.colorbar(shrink = 0.5)
plt.title("DEM Aspect", fontdict = font)
plt.show()
In [8]:
def hillshade(array, azimuth, angle_altitude):
        
    x, y = gradient(array)
    slope = pi / 2. - arctan(sqrt(x*x + y*y))
    aspect = arctan2(-x, y)
    azimuthrad = azimuth * pi / 180.
    altituderad = angle_altitude * pi / 180.
     
 
    shaded = sin(altituderad) * sin(slope)\
     + cos(altituderad) * cos(slope)\
     * cos(azimuthrad - aspect)
    return 255*(shaded + 1)/2
In [9]:
hillshade = hillshade(dem, 315, 45)
plt.imshow(hillshade, cmap = 'Greys')
plt.title("DEM Hillshade", fontdict = font)
plt.show()

Conclusion:


The Woosley Fire has caused many damages in the Los Angeles and Ventura Counties. Although the more urban and city areas were not ruined too much, they were in a very close proximity of the fire. This fire forced people to leave their homes and it was only a matter of time that the fire was contained. We can say that the vegetation will take a long time to be regained or it could possibly be permanently damaged. California state is already a very dry land so they already have high chances for wildfires to occur. Something to consider is water index for pre and post fire shows a very dramatic change. The state of California can use this finding to predict future fires because there is a pattern between the water index pre and post fire.

Future Plans / Implementations:

We believe from our findings, states like California can use water index and burn ratios to predict future and/or potential wildfires. They can also help create boundaries to protect homes and urban areas that are in high danger zones of fire. Or possibly have an evacuation plan ready for residents who are these potential wildfire zones. Another thing to consider as well is looking at land and finding ways to improve the land itself. That can prevent wildfires from happening and improve the ecosystem as well.

thank you for listening!

question / concerns?